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ABSTRACT 


Electromagnetic scattering from rough surfaces is of prime 
importance in the engineering field since it affects 
communications, radar, remote sensing, acoustics, etc. The 
actual problem of scattering from rough surfaces is 
complicated and involves three dimensional scattering from 
either lossy or dielectric, electrically large surfaces. 
Integral equations are widely utilized to solve this kind of 
problem but this solution to the problem is’ generally 
computationally intensive. 

In the On-Surface Radiation Boundary Condition (OSRC) 
method, a higher order radiation condition is imposed directly 
on the surface of the scatterer. This reduces the integral 
equation for the scattered field co a line integral which can 
be easily evaluated numerically. 

In this thesis, the OSRC method is used to formulate the 
problem of scattering from periodic rough, two-dimensional 
surfaces illuminated by a transverse magnetic, p'ane 
electromagnetic wave. Three geometric surfaces are 
considered. A comparison is made between the present 
formulation the exact solution, and the physical optics 


approximation. 
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I. INTRODUCTION 


A. HISTORY 

The problem of scattering from periodic rough surfaces has 
aroused the interest of physicists and engineers for the last 
70 years and has become of special interest in recent years, 
particularly in connection with the propagation of radio waves 
at frequencies above 30 MHz. The problem of scattering from 
rough surfaces arises, among others, in line-of-sight as well 
as mobile radio communication, radar, sea scattering, remote 
sensing, tropospheric propagation, radio astronomy, and 
acoustics. 

Mathematically, the problem of the scattering of waves by 
a rough surface is difficult to formulate and even more 


difficult to obtain numerical results. 


B. ON-8SURFACE RADIATION CONDITION (OSRC) METHOD 

Recently, a new formulation of electromagnetic scattering 
has been introduced ([Ref. 1). This formulation of 
electromagnetic wave scattering by convex, two-dimensional 
conducting bodies, the On-Surface Radiation Condition (OSRC) 
approach, is based upon an expansion of the radiation 
condition applied directly on the surface of a scatterer. 
According to the OSRC approach, the application of a suitable 
radiation condition directly on the conducting scattering 
surface can lead to substantial simplification of the 


frequency-domain integral equation for the scattered field. 
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For the transverse magnetic (TM) case, the integrand is known 


explicitly. For the transverse electric (TE) case, the 
integrand can be constructed by solving an ordinary 
differential equation around the scatterer surface contour. 

Implicit in the OSRC formulation is the fact that only 
outgoing waves are assumed to exist on the surface of the 
scatterer. Thus, the method will have serious drawbacks in 


cases where the latter condition is not met. 


C. STATEMENT OF THE PROBLEM 

In this thesis we apply the OSRC technique to scattering 
by a periodic rough’ surface. We consider a plane 
electromagnetic wave illuminating a two-dimensional, perfectly 
conducting periodic rough surface for the transverse magnetic 
(T™) polarization case (magnetic field parallel to the 
propagation plane or electric field parallel to the surface 
ridges as in Figure 1). 

The analysis of this problem includes the calculation of 
the induced electric current as well as the echo width of the 
scattering surfaces. Both of these quantities are calculated 
in Chapter II. 

In Chapter III, we calculate the current distribution and 
the echo width of a planar thin conducting strip, while in 
Chapter IV we calculate the same two quantities for a 
sinusoidal scattering surface and the surface current 


distribution on a scatterer of triangular shape. 











The numerical results of the current distributions as well 
as the echo width of the planar strip and the sinusoidal 
surface are presented in Chapter V. Comparisons of these 
results with results from the physical optics approximation 
and with the exact solution (method of moments), are also 
presented in this chapter. 

In Chapter VI, conclusions and recommendations are 
presented. Some mathematical derivations and the computer 
codes developed to evaluate and test our results are listed in 


the appendices. 





Ii. FORMULATION 


In this thesis we apply the OSRC technique to scattering 
by a rough surface. In this chapter we will formulate the 
problem and calculate both the surface current and the far- 
field echo width of a periodic rough, perfectly conducting 


surface. 


A. NOMENCLATURE 

We will consider a plane electromagnetic wave illuminating 
a two-dimensional, perfectly conducting, rough surface for the 
transverse magnetic (TM) polarization case (Figure 1). The 
width of the rough surface is considered finite and edge 
effects are ignored. 

We shall use Cartesian coordinates x, y, z with origin 0 
and unit vectors &, ¥, 2. The two dimensional rough surface 
is specified by the function z = f(x) where the mean level of 
this surface is the plane z = 0. 

All quantities, except angles, associated with the 
incident field will be dencted by the superscript i and those 


associated with the sca.tered field by the superscript s. 
Thus, the incident field is F and the scattered field £°. 
The medium above the surface is assumed to be free space. We 


shall assume E? to be linearly polarized. This implies that 


we need only the scalar value U* of the vector £?. 








B. GENERAL THEORY 


For TM polarization, the incident electric field £/ 


propagating perpendicular to the x-z plane at an angle 6, 
with respect to the +z axis, as shown in Figure 1, is given 


[Ref. 2] 
Bi = viesotp (2.1) 
where U? is the magnitude of the incident electric field 
Ui = @~Tho(x8in®, - 2008.) (2.2) 
ky is the wavenumber of the free space defined as 


20 anf, 
ee e 





and harmonic time dependence is assumed. The quantity A, is 
the free space wavelength associated with the incident 
electromagnetic wave of frequency f,, c is the speed of light 
in free space, and w = 2mf, is the radian frequency of the 


incident wave. 


The scattered electric field E*° is [Ref. 2] 
E*%=y%eietp (2.4) 


where U® is given, in the far-field, by [Ref. 1] 





i / 


-3(kor- =) Ss ° 
Ue = sl ou +jkycosbUt EFRON Ge! (2.5) 
[8k,nr]* © 
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Figure 1 Scattering geometry for a T polarized incident 
plane wave. 


where C represents the contour of the surface, A’ is the unit 
vector normal to the curve C at point P(x’,z’), as shown in 
Figure 2, and d@ /dn’ denotes an outward normal derivative on 
Cc. Note the integration in Equation (2.5) is with respect to 
the arc length, s', corresponding to the point P(x’,z’) at 
distance R from the origin, and that r is the distance of the 
observation point from the origin. The unit vector f in the 


direction of the observation point is 
£ = sinO- + cos6-2 (2.6) 


where 8 is the scattering angle, measured in opposite sense to 


the incidence angle @). 
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Figure 2 Notation for various vectors and angles. 


In Equation (2.5), the terms cosé and Rcosy are given by 


cosé = £:A/ (2.7) 

Reosy = R? (2.8) 
where 

R= x’'2+ 2/2 (2.9) 


The unknown quantity dU°/dn’ can be expressed in terms of 
the incident field using the second order radiation boundary 


condition. It is given by [Ref. 1] 


og 21 Ge): ik, 0? (x’) yie 1 PYF 02. 40) 
‘an’ | 2) 7°" BE -Fk,-C (NI 2[-7k)-C(x’)] as? 
where C(x’) is the curvature of the two dimensional 


scattering surface at point P(x’,z’) and a? jas’* is the 
second derivative with respect to the arc length of Cc. 
The curvature ((x) of the curve C that is defined by an 
equation of the type 2 = f(x) is given by [Ref. 3] 
d*z 
G(x) = dx" (2.11) 


a)" 


is shown in Appendix A to be 





while the quantity a*u' yas’? 








zu! -kjsin?6, + k3 2 cos8, (2sin®, - cose, 
= Ut 
ds” 1 +(2z\ 
= 


dz dz @z 
Fx 0380) Se dx dx 


b+ (Si 


By substituting Equations (2.11), (2.12), into Equation 


(sin@, - , 
2.12) 


(2.10) and carrying out some manipulation, we get 
































OUS _ yi dx? i- dx? + jk, 
an! 21 + dz\"]3 dz\7]4 . |dz 
ponodurst . 2 
| (a) | al jk,{t ates ax? 
1 
1+ dz a 2 
+ dx kj(sin?6, + cos6, Flesin®, - @ cos, \ 
dz\2)4 d? 
ofnfs » (Sey + Ss 
. dz daz 
(sine, - x0 ax 
+ jk, <2 cos®, + ax dx (2-13) 


The quantity 0U°/dn’ is expressed in terms of U/ and hence is 
completely known at all points on the rough surface. 

The unit vector fi’, normal to the curve C at point 
P(x’,z’), is given by [Ref. 3] 


dz 
f=ea/- 2% 2, 1» (2.14) 


1 1 
dz\?|45 dz\?|3 
1+ (== 1+{== 
where the upper sign corresponds to the outward normal and the 
lower sign to the inward normal. 
The formula for the partial derivative of the incident 


electric field with respect to the outward normal vector, 


which we will use later, is given by 














G2 sine, + cos8, 


au? = (Vu?) A! = jkui dx 
dn! 4s dz\2]+ 
lb + (ae) | 


(2.15) 


Substituting Equations (2.6) and (2.14) into Equation 


(2.7), we get 


~ 42 55n6 + cos8 


cos} = £1 = —————____—_— (2.16) 
ai 
dx 
Likewise substituting Equations (2.9) and (2.6) into 
Equation (2.8), we get 
k,Rcosy = k,R:£) = k,(xsin® + zcos6) (2.17) 


We will normalize all dimensions to the free space 


wavenumber k,. Accordingly, we let 





p=kx, t= kyz (2.18) 
so that 
os gh cae 
dx "%, 7 , az ao (2.19) 
and 
dz_ ada dz _, dt 
ax a ' aa * dp? eee 


Combining Equations (2.18), (2.19) and (2.20) with 
Equations (2.2), (2.16), (2.17) and (2.13), we express all 


quantities in terms of the normalized variables p,t as 


10 








yi = e J(psine, - tco88,) (2.2 1) 


--Z sino + cos8 
cos8 = (2.22) 


b+ (ST 


k,Rcoswy = psin®@ + tcos6é (2.23) 





at Sl 
2 
au* _ x yt d — dp? a4 


(SF) abi = psy es) 


dp? 
2 
1+ a : 2 
+ _ Petey sin*6, + $F cos8,{2sin8, - Si c0s8,} 














3 
: dt\*]5 , | d@t 
aljja + (SE) ]? + |S 
Fi (3s) | dp? 
sind, - = cos8, & 
+Jj a cos®, + a (2.24) 
aay, 


The arc length s is related to the surface z(x) by [Ref. 


3] 


11 








a 
ss fe + By)? ax (2.25) 
a 


so that the differential arc length ds’ can be expressed in 


terms of the normalized coordinates as 
! 1 dt\7]3 4-7 
a +: [2 2.26 
ds [1 | dp ( ) 


Substituting Equation (2.26) in Equation (2.5), we get: 





-I(ker- =) s jkpRcosy 274 

as e dU ° ; e at 2 / 2. 

ge || i + jkqcosb0'| a [2 + (SE) Jae (2.27) 
[8k xr]? © 


To organize our work better and to enable more efficient 


structural programming, we rewrite the above equation as 





U® = V(r) [F(p) -2(e) yap (2.28) 
¢ 
where 
~i{kor - 3) 
e 
V(r) = 7 (2.29) 
[8koRz] ? 
F(p) = + jk,cos8u! (2.30) 
12 





jk,Roose 2 
Op) = 2 _—_ [2 + Gaal (2.31) 


and 0U*/dn’ , cosé, u' and k,Rcosy as given by Equations 


(2.24), (2.22), (2.21) and (2.23) respectively. 


C. SURFACE CURRENT AND ECHO WIDTH 

Two quantities needed to compare our results with those in 
the literature are the electric surface current induced on the 
scattering surface and the far field echo width of this 
surface. These quantities will be presented next. 

Since vu'* = u'+ u* and the total magnetic field is given 


by 








Tot = 1 auret 
= 2.3 
psf Y= a ( 2) 


the normalized electric current on the scattering surface, 
equal to the ratio of the total to the incident magnetic 


field, is 





Hitot _ oben Oe" F au) (gi233) 


Hi jop,Hil on dn 


In the Equation (2.33) the only previously undefined term 
is the magnitude of the incident magnetic field, H', which can 
be obtained in terms of the electric field using the Maxwell's 


equation [Ref. 4) 
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For the case of a plane incident wave 


ky 


3 1 
eee ne 
0 (1) 





(2.35) 


where n, is the intrinsic impedance of the free space. 


The normalized surface electric current is now given by 


ret 1{oau* . gui P 1 {| oU% oui 
ai = 2{€ + au") = 304 $ + ar) (2.36) 








The far-field echo width o of the scattering surface is 


defined as [Ref. 2] 
o = lim, fe 2x1] (2.37) 


Since all quantities have been previously defined, o can be 


easily calculated. 
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III. THIN CONDUCTING STRIP 


The mathematical derivations of the previous chapter are 
for an arbitrary rough surface. In this chapter we specialize 
to the case of TM scattering by a thin conducting strip. 
Because results for the induced current and echo width for a 
conducting strip are known by other means this will serve as 


a check to our formulation. 


A. SURFACE CURRENT 
For a thin planar perfectly conducting strip with the 


geometry of Figure 3, the surface is described by 
Z=ax , t=ap (3.1) 


where a is the slope of the surface. 

Using the results of the previous chapter, the normalized 
surface electric current density on this surface is given by 
( _ 8in?8, +acos6,(2sin8, eee) , 25in6, +cos8, 


Tot 
A” sul 
H? 2(1 + a?) 2 


| (3.2) 
[. +a]? 





Notice that the current distribution is known analytically. 
This is in contrast to other methods where differential or 
integral equations must be solved to obtain the current 
distribution. 

Let us now derive a formula for the scattered echo width. 


According to Appendix B, we can express the integral in 


15 











Figure 3 Planar, thin conducting strip 
illuminated by a T™ polarized, plane wave. 


Equation (2.28) as 


pa+kol 


[F(p)-o(p)dp = 2 f oe O(p) dp 
¢ pe-Kob 


Equation (2.27) becomes 


a¢ 





k,u4 on 


profile, 


(3.3) 


Kol 
ute =vir)2 f (z z ia | sees +(£) Pap (3.4) 


pe-Kol 


Note that no special consideration for the fields is given at 


the edges x = + L. 
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Other quantities needed to evaluate the integral are 








yi = oe Je(eine, - acos6,) (3.5) 
k,Rcosy = p(sin6 + acos6) (3.6) 
1 ous 1 
=jil - ———D 3.7 
kU! on { 2] + a7] mee 
where 
D = sin*6, + acos6,(2sin®, - acos®,) (3.8) 


Substituting the above equations into Equation (3.4), we 


get 
+koL 
U* = V(r) 2ifa - aah + a’? f eJ* dp (3.9) 
2[1 + a7] koh 
where 
G = (sin® - sin6,) + a(cos®8 + cos6,) (3.10) 


The integral in Equation (3.9) can be evaluated in closed 
form. Carrying out the integration and substituting for V(r), 


we get 





of 2 4a? 4 ees D__), | 8im(koLG) -i(ker-=) 3.22) 
: 2[1 + a?]} ° : 


8k, "Ir k,LG 
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B. FAR-FIELD ECHO WIDTH 

Now that the scattered field is known, the echo width can 
be easily calculated using Equation (2.37). Substituting 
Equations (2.2) and (3.11) into Equation (2.37), and 
normalizing the echo width o to the free space wavelength Aj, 


we obtain 


(3.12) 





2 
; : 
go tates - 2a) ae 

The above formula is for the bistatic echo width (receiver 

and transmitter are considered at different locations). 
The monostatic case (transmitter and receiver at the same 
location) is a special case of the bistatic case, where the 
angle of incidence 8, is related with the observation angle 


8 by 


6 = 2x - 0, (3.13) 
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IV. SCATTERING FROM TRIANGULAR AND SINUSOIDAL 
PERFECTLY CONDUCTING SURFACE 


In this chapter we investigate TM scattering from 
perfectly conducting triangular and sinusoidal surfaces. 
Expressions for the electric current induced for both the 
surfaces and the echo width for the sinusoidal case alone are 
developed. The sinusoidal surface is generally considered as 
a prototype surface for studying scattering from more 
complicated, rough surfaces. The scattering surface profiles 


are shown in Figure 4. 


A. TRIANGULAR CASE 
For the case of the triangular profile, as shown in Figure 


4(a), the surface is described by 


: _ 4h d 
z(x) = A q* ; Osx< 3 (4.1a) 
4h d 
z(x) = -3h + “x 5 <xsd (4.1b) 


where d is the wavelength of the periodical surface. The 
expression for the normalized surface electric current density 


(Equation (2.36)) is simplified to 
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Figure 4 Profiles of periodic surfaces: (a) Triangular; 
(b) Sinusoidal. 
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Tet ae sin’, +  cos®,{2sin8, - 2 cos8, 
H? : ; 2j1 + (2) 
ax 
Z sinb, + cos8, 
+ —_ 
1 + (Sy 3 a 
dx 
where 
@-- ’ osxs¥ (4.3a) 
daz 4h d 
we og. it Se SS (4.3b) 


B. SINUSOIDAL CASE 
The sinusoidal scattering surface profile (Fig. 4(b)) is 


described by 


z(x) = hcos (wx) (4.4) 
where w = 2H7/d and da is the wavelength of the sinusoidal 
surface. Hence 

az . -hwsin(wx) , az . -hw*cos (wx) (4.5) 








and 


k,z = k,hcos(wx) = +t = kghcos{ #2) (4.7) 
0 


These quantities are used to express the normalized 
surface current induced on the surface and defined in Equation 
(2.36). The echo width of this surface is obtained in the 
following fashion. 

Since the scattering surface is perfectly conducting with 
considerable thickness, the integral in Equation (2.28) is 


shown in Appendix B to be 


prt+kyoL 
ou* . gui 
F(p):O(p) dp = (az + a") 9(p) dp (4.8) 
i) ay on on 


The far-zone, scattered field reduces to 





P=*Ko aus aut 
ore vies f [aaa « 2) 


: i seaaeas ; (S£)"]Fap (4-9) 
penkytl Ko? on dp 





where 


k,Rcosy = psin® + kghcos| #2 Jcose (4.10) 


0 


The incident field on the surface of the scatterer is 


Hive Bs ase - kyhcos{ FP conte) (4.11) 


For computational purposes it is convenient to rewrite the 


incident field as 
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yigikoRcost _ og IT (4.12) 


where 


T = p(siné - sin6,) + kghcod $2 | (cost - cos6,) (4.13) 
Q 


Substituting these in Equation (4.9), we obtain 


prekyL 
u* = vir) f H(p) dp (4.14) 
pe-Kol 
where 
= —1_[9U* , 9U*),5 dt \7)3 
se sal an ~ On Je p “lee eapae? 


We normalize the limits of integration to the quantity k,-L. 


After some manipulation, Equation (4.14) reduces to 


Uf=V(r) 








0 * Kol . . 
[ H(prdp + f H(p) dp] = oF = viz) [vt(p) +H(-p)ldp (4.16) 
~keL 0 0 


With the substitution p = k,Lu , we get 


U* = V(r) koLM (4.17) 


where 
1 
M = f[H(k.Lu) + H(-K,Lu) jdu (4.18) 
i] 


Since a closed form formula has not been found for the 
scattered field, the integration in Equation (4.18) is 


performed numerically. 
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Once the quantity M has been calculated, the echo width o 


normalized to the wavelength A, is obtained from 


Oo «| tye (Koh)? 
ie |M| a (4.19) 


When normalized with respect to the length (2L) of the 


scattering surface, the echo width is 


se = |m|?—2— (4.20) 


kL 
2L 8 
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Vv. NUMERICAL RESULTS 


In this chapter, computed results for the induced 
electric surface current and the echo width are presented. 
In order te check our theory and computer codes, the results 
of the analytical solution (OSRC) are compared against the 
moment method, treated here as exact, and the physical optics 
solution. Details of these comparisons are presented in the 
following sections. 

In all models where numerical results are presented, the 
following factors are considered: 

1. The incident electromagnetic wave is a Transverse 
Magnetic plane wave incident in the x-z plane. 

2. However, in an actual periodic surface the surface is 
infinite in the direction of periodicity, the width of the 
scattering surface considered here is finite and edge effects 
are ignored. 

3. The induced electric current is plotted against the 
normalized distance (across one period of the periodic 
surface). 

4. The resulting echo width o is normalized to the 
wavelength A, for the planar strip and to the surface total 
length 2L for the sinusoidal surface. 

5. Curves presented for the exact and physical optics 


solution are taken from [Ref. 5] and [Ref. 6]. 
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A. SURFACE CURRENT COMPARISONS 

For a planar strip (Figure 3) of length 2L, the magnitude 
of the surface current has a uniform distribution while the 
phase drops linearly across the strip, as shown in Figure 5, 
for incidence angles 6 = 0 ,30, 60, 90. Notice that the 
current distribution does not exhibit any special behavior at 
the edges. 

The magnitude of the surface current as a function of the 
incidence angle 8, for the OSRC solution as compared with the 
physical optics solution is shown in Figure 6. We can see 
that the agreement is remarkable between the two techniques 
for angles of incidence close to zero. However, the results 
deviate for angles of incidence exceeding 45°. For normal 
incidence the numerical (OSRC) results agreed with the exact 
solution to an accuracy better than 0.01 percent. 

The comparisons for the sinusoidal surface between the 
numerical results generated by the present approach with both 
the exact solution and the physical optics approximation are 
shown in Figures 7, 8 and 9. Data pertaining to the latter 
two techniques are taken from [Ref. 5]. 

The surface currents obtained by all three methods for 
normal incidence are shown in Figures 7 and 8. Agreement of 
our solution with the other two is best for h = 0.25A, and 
d = 3.94, (Figure 8), while for h = 0.254, and d = 1.9A, 


(Figure 7) the numerical results deviate considerably from the 
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other two. The results seem to be closer to the physical 
optics approximation in both cases. The maximum deviation 
occurs at the peak and trough of the surface where the 
curvature is maximum. This is expected of the OSRC technique. 

In general both outgoing and incoming waves exist on the 
surface of the scatterer. Millar {Ref. 7} has shown that for 
2mh/d < 0.448, the fields on the surface are predominantly of 
outgoing type. In the OSRC technique, a higher order 
radiation boundary condition is imposed directly on the 
surface of the scatterer. This implies that only outgoing 
waves have been assumed to exist on the surface of the 
scatterer. In the case of Figure 8, 2mh/d = 0.403, which is 
less than the critical value of 0.448. Hence better agreement 
is seen for this case. However, in Figure 7, 2mh/d = 0.826 
resulting in poorer agreement. 

In Figure 9 we present the results for the same geometry 
as in Figure 8, but for a 75° angle of incidence. The OSRC 
solution is not affected by the occurrence of shadowing as in 
the case of the physical optics. However, the results differ 
from the exact results especially in the central region 
(shadow region). 

The current distributions on periodic surfaces with 
sinusoidal and triangular profiles (see Figure 4) are shown in 
Figures 10 through 15 represent. Each surface has a period 


d = 0.2A,, and a maximum amplitude h = 0.14, (2mh/d =m). 
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The angles of incidence considered are 0) = 0°, 30° and 60°. 


In each figure, the OSRC solution results are shown along with 
the exact solution as given in (Ref. 6]. It is seen from 
Figure 10 that the OSRC technique predicts an amplitude of 
approximately 6, both at the crest and the trough of a 
sinusoidal profile, whereas the exact solution predicts values 
of 5 and 0 at the crest and trough, respectively. The rather 
large deviation of the results at the trough is attributed to 
the presence of standing waves, which the OSRC technique fails 
to take into account. 

For oblique incidences, OSRC technique predicts a jump 
discontinuity in the current magnitude at the center for 
triangular profiles as observed in Figures 13 through 15. 
This is due to the discontinuity of the surface slopes at 
these locations. The exact solution exhibits the proper edge 
condition whereas the OSRC fails to do so. A similar jump 


discontinuity is predicted by the physical optics [Ref. 6]. 


B. ECHO WIDTH COMPUTATIONS - COMPARISONS 

The normalized monostatic echo width o/A, of a plane strip 
with length 2L, as a function of the incidence angle 8, is 
shown in Figure 16. The normalized bistatic echo width o/d), 
as a function of the observation angle 6, for kL = 5, kL = 
10 and angle of incidence 0° and 60° is shown in Figure 17. 

Figure 16 shows also the comparison between the numerical, 
the exact (Ref. 1], and the physical optics [Ref. 8) results. 
We see that we have a very good agreement (within 1 dB) with 
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the physical optics solution for angle of incidence up to 70°. 
Agreement with the exact solution is poor beyond 30° angle of 
incidence. The deviation from the exact solution, for larger 
angles of incidence, is probably due to edge currents that 
have been ignored in the formulation. 

The normalized monostatic echo width o/2L, of a sinusoidal 
surface of length 2L with parameters h = 0.254, and d = 3.9A, 
is shown in Figure 18 while Figure 19 shows the normalized 
bistatic echo width o/2L of a sinusoidal surface with the same 
profile, for angle of incidence 0° and 60°. No comparison is 


made for the echo width of the sinusoidal surface. 
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(b) 
Figure 5 Normalized surface current versus distance qcross 
a horizontal strip, L = 0.24,. Lines: (——) - 6) = 9; 
(-+-) - 6, = 30; ( x-) ~- 6, = 60; (:o-) - 6, = 90. (a) 
Magnitude; (b) Phase. 
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Figure 6 Normalized surface current versus angle of 
incidence 6, across a horizontal thin strip. Solid line: 
OSRC results; broken line: physical optics approximation. 
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(b) 
Figure 7 Normalized surface current, versus distance across 


period of sinusoidal surface, 6) = 0, d = 1.94), h = 0.25A). 
Solid lines: OSRC results; broken lines: (-x-) - exact, 


(.-o- :) - physical optics. 
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Figure 8 Normalized surface current, versus distance across 
period of sinusoidal surface, 8, = 0, d = 3.9A,, h = 0.25A). 
Solid lines: OSRC results; broken lines: (-x-) - exact, 

( ~O- :) - physical optics. 
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(b) 
Figure 9 Normalized surface current yersus distance across 
period of sinusoidal surface, 6, = 75, d= 3.94,, h = 
0.25A4,. Solid lines: OSRC results; broken lines: (-x-) - 
exact, (-~-0--:) - physical optics. 
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(b) 
Figure 10 Normalized surface current versus distance 


across period of sinusoidal surface, 8, = 0, d = 0.2A), 
h = 0.1A,. Solid lines: OSRC results; broken lines: exact 


solution. 
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Figure 11 Normalized surface current versus, distance 


across period of sinusoidal surface, 6, = 30, d = 0.2A,, 


h = 0.1A,. Solid lines: OSRC results; broken lines: 
solution. 
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Figure 12 Normalized surface current versus distance 
across period of sinusoidal surface, 0, = 60, d = 0.2d), 
h = 0.1A4,. Solid lines: OSRC results; broken lines: exact 


solution. 
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(b) 
Figure 13 Normalized surface current versus distance 
across period of triangular surface, 6, = 0, d = 0.2A), 


solution. 
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Figure 14 Normalized surface current versu 


across period of triangular surface, 6, = 30, d= 0.2A5, 
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h = 0.14,. Solid lines: OSRC results; broken lines: 


solution. 
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Figure 15 Normalized surface current versu 
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Solid lines: OSRC results; broken lines: 
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Figure 16 Normalized monostatic echo width o/a,, of 
horizontal conducting strip. Solid lines: OSRC’ results; 
broken lines: (+«x-) - exact solution, (---) - physical 
optics. (a) kL = 5; (b) kL = 10. 
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Figure 17 Normalized bistatic echo width o/A, of 
horizontal conducting strip - Numerical analysis (OSRC) 
results. Solid lines: kL = 10; broken lines: kL = 5. 


(a) 8,=0°; (b) @,=60°. 
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Figure 18 Normalized monostatic echo width o/2L of 
sinusoidal surface, d = 3.94), h = 0.2545; - OSRC results -. 
Solid line: kL = 10; broken line: kL = 5. 
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Figure 19 Normalized bistatic echo width o/2L of 
sinusoidal surface, d = 3.9A,, h = 0.254); - OSRC results ie 
Solid lines: kL = 10; broken lines: kL = 5. (a) 0@,= 0; 
(b) 6, = 60. 
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VI. CONCLUSIONS 


A. SUMMARY 

Electromagnetic scattering from rough surfaces is hard to 
formulate and even harder to solve. In the On-Surface 
Radiation Boundary Condition a higher order radiation 
condition is imposed directly on the surface of the scatterer. 
The frequency-domain integral equation for the scattered field 
is reduced to just a line integral. This leads to substantial 
simplification of the problem. Analytical expressions are 
obtained for the surface scattered fields. 

In this thesis, the OSRC method has been used to formulate 
electromagnetic scattering from a periodic rough, perfectly 
conducting surface. Transverse magnetic (TM) polarized, plane 
incident electromagnetic wave was considered throughout the 
whole formulation; the surface was considered of finite 
length, and edge effects were ignored. 

Three different scattering surface profiles have been 
considered to test our formulation. The surface currents as 
well as the echo width of these surfaces were calculated by 
developing computer codes and the results were compared with 
results from the exact solution and the physical optics 
approximation. 

For all the examined surfaces, the OSRC method predicts 


results somewhat between the two other solutions. For the 
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case of a thin strip, the OSRC results for the echo width are 
within 0.1 dB (except at nulls) of the exact solution for 
angles of incidence between o° and 30° and deteriorating for 
larger angles of incidence. Actually OSRC results are more 
similar to the results given by physical optics. 

For the cases of a triangular surface and a sinusoidal 
surface with parameters h and d such that 2mh/d > 0.448, the 
OSRC results are not very accurate. 

For the case of a sinusoidal surface, the validity of the 
OSRC results was shown to depend on surface profile parameters 
as well as the angle of incidence. When the quantity 2m7h/d 
< 0.448 and the angle of incidence is close to zero, the OSRC 
results are very accurate, while for larger angles of 
incidence it is somewhat closer to the exact solution than the 
physical optics one. Overall, the results of the OSRC 
technique seem to mimic those of the physical optics rather 
than the exact solution. This is clearly seen in all 


comparisons of current distribution and echo width. 


B. RECOMMENDATIONS 

With the basic theory for the TM case tested, we recommend 
the investigation of surface currents and echo width of other 
periodic surfaces. The formulation and the validation of the 
same problem, for the Transverse Electric (TE) case is another 
recommendation. Finally, how much the edge and the finite 


surface length affect the results needs to be investigated. 
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APPENDIX A 
MATHEMATIC DERIVATIONS 
In this Appendix we will derive the second partial 
derivative of the incident field with respect to the arc | 
length on the scattering surface. 
Let the electric field of the incident wave be denoted by 
an arbitrary function of the type U = U(x,z). The surface cf 


the scatterer is described by z = z(x). 


au _ aUdx , Adz _/aU , aU dz) dx 
ds oxds dz ds (5. * Sz ae nen) 


where s is the arc length on the scattering surface. 


Since 
ds _ dz\")3 A.2 
ae [2 + (¥)] ( ) 
we have then 
oU QU dz 


dU _ Ox dz dx 
Sep ce ee (A.3) 
Zz 





See Ae oa ee (A. 4) 
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It can be shown that 





fu, PU dz , Ww dz oU , OU dz\ dz d’*z 
ow _ ax? * Oxdz “ax” OZ dx? tga Oz dx} dx dx? (A.5) 


“bleh (ey 


and 
eu , Pu dz 
2 
Le OZzOx Oz = (A.6) 
© b+ (aT 
dx 
For the case of a plane incident wave 
is e Io (xsin®y-z00S8,) (A.7) 
we have 
r? : Z i ; i 
oe = -jk,e I(kx - kz) _ -jk,ui (A.8) 
OU". ~jk,e Ut -**) = jk,U (A.9) 
Gs ee = JX; : 








i . 

a = -k2e Whe - Ket) - 42 yi (A.10) 
i : 

— = k2e Ihe - ket) 2 2 yi (A.22) 
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@ui ui = ~j(k,X - k,2) 
Oxdz ~ azdx ~ *xke° ore 





where 
k, = k,sin6, , k, = k,cos8, (A.13) 


Finally, substituting Equations (A.5) through (A.13) in 


Equation (A.4), we get 








-kgsin?6,+ké 22 oos8 2sin6, - 2 cos6, + jk, 2-088, 
Pui _ yt °'dx dx dx? 
ds? dz\? 
1+ (= 2 
: dz dz d*z 
(sin@, - —cos8,) = 
+ jk, dee ax dx? (A214) 


tT 
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APPENDIX B 
EVALUATION OF INTEGRALS 

In this Appendix we present an approach for evaluating the 
line integrals with respect to the arc length that are encountered 
in Chapters 3 and 4 for the case of a planar strip and sinusoidal 
surface. 

We are interested in evaluating an integral of the form 
J-( )dp where C is a conformal contour enclosing the surface. 

We first consider the planar strip. The closed contour C 
around the scattering surface is considered to comprise of four 
contours C,, C,, C; and C,, as shown in Figure Bl. 

The contribution from contours C, and C, approaches zero for 
bounded integrands. The arc length s, as well as the differential 
arc length, ds, on the two contours C, and C,, with respect to the 


point A (Figure Bl) is 
_ -£ 
B= 5+xX , ds = dx (B.1) 
on C,, and 
= L 
s-L+(5-x) , ds = -dx (B.2) 


on C,. 
For computational purposes let us name the integral in 


Equation (2.5) as 





I= {= + jkuicoss }as (B.3) 
¢ 
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Figure Bl. Integration contour for a thin strip. 


We then have 
Ls eon dae * | ae Pal eaae 2 aae) (B.4) 


where fi, and fA, are the unit normal vectors outward to the contours 


C, and C, respectively. 


Using Equations (B.1) and (B.2), we see that Equation (B.4) 


becomes 
. pause us 
re [ax - fF ian ae a poe) (B.S) 
c C3 
Since Cc, = -C, and fi, = -A,, we have 
s dus 
r= flax - f 2 Bg IR | reas | eae) (B.6) 
Cy “cy “cy 
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7 fea + [ax , kyo fetax - ae (B.7) 
Cy Cy C, cy 


Hence, for the case of a thin conducting strip, we have 





{= + jkv!coss ds = 2f Zax (B.8) 
c Cy 


Let us now examine the case where the thickness of the 
scattering surface is no longer small, so that the second contour 
C, is considered to lie inside the conducting surface (Figure B2). 


In this case we have 


ous 
Is mn ——ds + jk, on f aid (B.9) 


Cy +g Cy +C2 


The second term in the above Equation (B.9) vanishes on the closed 
contour C,+C, as shown previously. Inside a perfect conductor, 


the total magnetic field is zero. Hence on the contour C,, ou’/dan 


= -du'/an. Substituting these, we see that Equation (B.9) 
becomes 
ou* oui 
r= [5 -4s + PE J an, (B.10) 
¢ 2 Cc, 2 
Since C, = -C, and fi, = -fA, = fi, it follows that 

















ous 1 Ce} 8 fad 
ff an +jk,U cost |ds = Jan, 





_ f{ ous , ou? 
ds -{{ ant so a0" (B.11) 
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Figure B2. Integration contour for a thick surface. 
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APPENDIX C 
PROGRAMS 
In this Appendix we present the computer programs, all in 
FORTRAN language that are needed to perform the various 


calculations. 


A. DISCUSSION 

Programs STRIPCUR, SINUSCUR and TRIANCUR calculate the 
surface currents induced on the scattering surfaces, for a 
strip, sinusoidal and triangular profiles, respectively. 

Programs STRIPMON and STRIPBI calculate the normalized 
monostatic and bistatic echo width, respectively, of thin, 
planar, perfectly conducting strip. Programs SINUSMON and 
SINUSBI calculate the same quantities but for a sinusoidal 
surface. 

Since we deal with complex numbers most of the above 
mentioned programs contain subroutines to calculate the real 
(i.e. STRPREAL, SINREAL, FRSTREAL etc.) and the imaginary part 
(i.e. STRPIMAG, SINIMAG, MOODIMAG etc.) of a quantity and use 
subroutine PHASER to calculate the phase of this complex 
quantity. 

Program TRIANCUR calculates the current distribution on 
the triangular surface, by examining separatly the two halves 
of the surface, the first with negative slope (subroutines 
FRSTREAL, FRSTIMAG), and the second with positive slope 
(subroutines SCNDREAL, SCNDIMAG). 

Programs SINUSMON and SINUSBI carry out the integration of 
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Equation (4.18), using Filon's method, as it is given by 
function subprogram FILON [Ref. 9]. To take advantage of the 
fast and acurate results of function FILON, we had to separate 
the even and the odd part of the examined function, as they 
are shown in subroutines MOEVREAL, MOEVIMAG, BIEVREAL, 
BIEVIMAG, and MOODREAL, MOODIMAG, BIODREAL, BIODIMAG 


respectively. 
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gaaNAANANANANNANNANNNNNNANNANANNANNANANANANANnNAN 


PROGRAM STRIPCUR 


THIS PROGRAM CALCULATES THE NORMALIZED SURFACE CURRENT 
ON A PLANAR, PERFECTLY CONDUCTING STRIP, DESCRIBED BY 
THE EQUATION: z(x) =A * x 

ATA: 
F = FREQUENCY IN GHz 


INCID = ELECTRIC FIELD ANGLE OF INCIDENCE (IN DEGREES), 
WITH RESPECT TO THE VERTICAL (+Z) AXIS 


A = SLOPE OF THE STRIP 


= CONSIDERED LENGTH OF THE SURFACE 


M = NUMBER OF EQUIDISTANT POINTS BETWEEN x=0 AND 
x=L WHERE OUR FUNCTION IS GOING TO BE 
EVALUATED. 

OUTPUT DATA: 


FREAL,FIMAG = REAL AND IMAGINARY PART RESPECTIVELY OF 
THE NORMALIZED SURFACE CURRENT 


AMAGNIT, PHASE = MAGNITUDE AND PHASE RESPECTIVELY OF 
THE NORMALIZED SURFACE CURRENT 


REAL K, INCID, INCR, L 

COMPLEX ZET 

COMMON K, A, THETAO 

OPEN (UNIT=11, FILE='STRIPCUR.DAT', FORM='FORMATTED', 
STATUS=' UNKNOWN ' ) 

READ(5,*)  F,INCID,A,L,M 


PI = 4.*ATAN(1.) 

RAD = PI/180. 

THETAO = INCID*RAD 

K = 20.*PI*F/3. 

INCR = L/M 

DO 100 X = 0., L, INCR 
DISTANCE = X/L 
FREAL = STRPREAL(X) 
FIMAG = STRPIMAG(X) 
ZET = CMPLX(FREAL, FIMAG) 
AMAGNIT = CABS (ZET) 
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PHASE = PHASER(FREAL, FIMAG) 
WRITE(11,20) DISTANCE, FREAL, FIMAG, AMAGNIT, PHASE 
FORMAT (1X, F12.6, 4X, F12.6, 4X, F12.6, 4X, F12.6, 
4X, F12.6) 
CONTINUE 
STOP 
END 


REKKKEKKRKEKKKERKEKKKRRKKKKKKES 


REAL FUNCTION STRPREAL(X) 
REREKEKHREKEREKREEREKERKKEKE 


REAL K, KC3 
COMPLEX J, POWER, FUSCAT, FUINC, UI, FUNC 
COMMON K, A, THETAO 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 


Z = A*X 
J = CMPLX(0.,1.) 
POWER = - J*K*(X*STHO - Z*CTHO) 


C = SQRT(1.+ A*A) 

D = (STHO)**2.+ A*CTHO*(2.*STHO - A*CTHO) 
FUSCAT = J*(1. - D/(2.*C#C)) 

FUINC = J*((A*STHO + CTHO)/C) 

UI = CEXP(POWER) 

FUNC = - J#UI*(FUSCAT + FUINC) 

FUNCREAL = REAL(FUNC) 

RETURN 

END 


RREKKKKKKKRKKKKKKKKKKRKRKKKE 


REAL FUNCTION STRPIMAG(X) 
RRKEAEKRRERKREKEEEEEEREREREEE 


REAL K, KC3 
COMPLEX J, POWER, FUSCAT, FUINC, UI, FUNC 
COMMON K, A, THETAO 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 

Z = AkX 

J = CMPLX(0.,1.) 

POWER = - J*K*(X*STHO - Z*CTHO) 

C = SQRT(1.+ A*A) 

D = (STHO)*#*2.+ A*CTHO*(2.*STHO - A*CTHO) 
FUSCAT = J#(1. - D/(2.*C#C)) 

FUINC = J*((A*STHO + CTHO)/C) 

UI = CEXP(POWER) 

FUNC = - J#UI*(FUSCAT + FUINC) 

FUNCIMAG = AIMAG(FUNC) 

RETURN 








Cc 
Cc 
Cc 
Cc 
Cc 








RAERKREKAEARRKKREEKREERKKKERREEEEK 


REAL FUNCTION PHASER(X, Z) 
RHEKKKEKEKKKKKKKKKKEKKKEEKKE 


PI=f 
PI = 4.*ATAN(1.) 


IF (X.EQ.0.0.AND.2.EQ.0.0) THEN 
ANGLE = 0. 
ELSEIF (X.EQ.0.0.AND.Z.GT.0.0) THEN 
ANGLE = PI/2. 
ELSEIF (X.EQ.0.0.AND.Z.LT.0.0) THEN 
ANGLE = ~-PI/2. 
ELSEIF (X.GT.0.0.AND.Z.GE.0.0) THEN 
ANGLE = ATAN(ABS(Z/X) ) 
ELSEIF (X.GT.0.0.AND.Z.LT.0.0) THEN 
ANGLE = - ATAN(ABS(Z/X)) 
ELSEIF (X.LT.0.0.AND.Z.LT.0.0) THEN 
ANGLE = -PI/2. - ATAN(ABS(Z/X) ) 
ELSEIF (X.LT.0.0.AND.Z.GE.0.0) THEN 
ANGLE = PI - ATAN(ABS(Z/X)) 





ENDIF 
COST = 180/x 

CONST = 57.29577951 

PHASER = ANGLE*CONST 

RETURN 

END 














Cc. 
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PROGRAM SINUSCUR 


THIS PROGRAM CALCULATES THE NORMALIZED SURFACE CURRENT 
ON A SINUSOIDAL, PERFECTLY CONDUCTING SURFACE, DESCRIBED 


BY THE EQUATION: z(x) = H*cos([2*pi/d]*x) 
NP ATA: 
F = FREQUENCY IN GHz 


INCID = ELECTRIC FIELD ANGLE OF INCIDENCE (IN DEGREES) 
WITH RESPECT TO THE VERTICAL (+2) AXIS 


H, a = PARAMETERS OF THE SINUSOIDAL SURFACE 

M = NUMBER OF EQUIDISTANT POINTS BETWEEN X=0 AND 
=d WHERE OUR FUNCTION IS GOING TO BE 
EVALUATED 

OUTPUT DATA: 


FREAL,FIMAG = REAL AND IMAGINARY PART RESPECTIVELY OF 
THE NORMALIZED SURFACE CURRENT 


AMAGNIT, PHASE = MAGNITUDE AND PHASE RESPECTIVELY OF 
THE NORMALIZED SURFACE CURRENT 


REAL K, INCID, INCR 
COPLEX ZET 
COMMON K, H, W, THETAO 


OPEN (UNIT=12, FILE='SINUSCUR.DAT', FORM='FORMATTED', 
STATUS=' UNKNOWN ' ) 
READ(5,*) F, INCID, H, a, M 


I 


P 4.*ATAN(1.) 
w 2.*PI/d 

RAD = PI/180. 

THETAO = INCID*RAD 

K = 20.*PI*F/3. 

INCR = a/M 

DO 100 X = 0., d, INCR 
DISTANCE = x/d 

FREAL = SINREAL(X) 

FIMAG = SINIMAG(X) 

ZET = CMPLX(FREAL, FIMAG) 
AMAGNIT = CABS(ZET) 

PHASE = PHASER(FREAL, FIMAG) 














100 


QaaQq AaNQ 


aq an 





WRITE(12,20) DISTANCE, FREAL, FIMAG, AMAGNIT, PHASE 
FORMAT(1X, F12.6, 4X, F12.6, 4X, F12.6, 4X, F12.6, 


4X, F12.6) 
CONTINUE 
STOP 
END 


KkkkkkKKKRRKKKKEKKEKKKEKKEE 


REAL FUNCTION SINREAL(X) 
REREKEKRREKKEEEKEKEKKKKKEE 


REAL K, KC3 
COMPLEX J, POWER, BJKC3, DJE, FUSCAT, FUINC, UI, 
COMMON K, H, W, THETAO 


STHO = SIN(THETAO) 
CTHO = COS (THETAO) 


WX = Wex 
SW = SIN(WX) 
CW = COS (WX) 
Z = HeCW 


J = CMPLX(0.,1.) 

POWER = - J*K*(X*STHO - Z*CTHO) 

A = -H*WesSw 

B = -H*WeW*CW 

BA = ABS(B) 

C = SQRT(1.+ A*A) 

D = K*((STHO) **2.+ A*CTHO*(2.*STHO - A*CTHO)) 

E = B*(CTHO + A*(STHO - A*CTHO)/(C*C) ) 

KC3 = K*C#C#C 

BJKC3 = BA + J*KC3 

DIE =D + J#E 

FUSCAT = BA/(2.*KC3)*(1 - BA/(4.*BJKC3)) + J 
+ C/(2.*BJKC3) *DJE 

FUINC = J*((A*STHO + CTHO)/C) 

UI = CEXP(POWER) 

FUNC = - J*#UI*(FUSCAT + FUINC) 

FUNCREAL = REAL(FUNC) 

RETURN 


RREAEARRRKARARKEREREERRRKRKKRRKE 


REAL FUNCTION SINIMAG(X) 
REAREKREREREEREEEEKREREREE 


REAL K, KC3 
COMPLEX J, POWER, BJKC3, DJE, FUSCAT, FUINC, UI, 
COMMON K, H, W, THETAO 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
WX = W*X 


FUNC 


FUNC 











SW = SIN(WX) 

CW = COS(WX) 

Z = H*CW 

J = CMPLX(0.,1.) 

POWER = - J*#K*(X*STHO - Z*CTHO) 
A = -H*W*sSw 

B = -H*W*Ww*Cw 

BA = ABS(B) 


C = SQRT(1.+ A*A) 

D = K*((STHO)**2.+ A*CTHO*(2.*STHO - A*CTHO) ) 

E = B*(CTHO + A*(STHO - A*CTHO)/(C#C) ) 

KC3. = K*Cecec 

BJKC3 = BA + J*KC3 

DIE =D + J*E 

FUSCAT = BA/(2.*KC3)*(1 - BA/(4.*BJKC3)) + J 
+ C/(2.*BJKC3) *DJE 

FUINC = J*((A*STHO + CTHO)/C) 

UI = CEXP(POWER) 


FUNC = - J*UI*(FUSCAT + FUINC) 
FUNCIMAG = AIMAG(FUNC) 
RETURN 


END 








D. 
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PROGRAM TRIANCUR 


THIS PROGRAM CALCULATES THE NORMALIZED SURFACE CURRENT 
ON A TRIANGULAR, PERFECTLY CONDUCTING SURFACE, DESCRIBED 
BY THE EQUATION: 


z(x) = H - (4*H/d)*x , for 0 < x < @/2 (FIRST HALF) 
z2(X) = -3*H + (4*H/d)*x , for d/2 < x < d (SECOND HALF) 
INPUT DATA: 

F = FREQUENCY IN GHz 


INCID = ELECTRIC FIELD ANGLE OF INCIDENCE (DEGREES), 
WITH RESPECT TO THE VERTICAL (+Z) AXIS 


H, a = PARAMETERS OF THE TRIANGULAR SURFACE 

M = NUMBER OF EQUIDISTANT POINTS BETWEEN x=0 AND 
x=d WHERE OUR FUNCTION IS GOING TO BE 
EVALUATED 

OUTPUT DATA: 


FREAL,FIMAG = REAL AND IMAGINARY PART RESPECTIVELY OF 
THE NORMALIZED SURFACE CURRENT 


AMAGNIT, PHASE = MAGNITUDE AND PHASE RESPECTIVELY OF 
THE NORMALIZED SURFACE CURRENT 


REAL K, INCID, INCR 

COMPLEX ZET 

COMMON K, H, d, THETAO 

OPEN (UNIT=14, FILE='TRIANCUR.DAT', FORM='FORMATTED', 
STATUS=' UNKNOWN ! ) 

READ (5, *) F, INCID, H, da, M 


PI = 4.*ATAN(1.) 
RAD = PI/180. 
THETAO = INCID*RAD 
K = 20.*PI*F/3. 
INCR = d/M 
DO 100 X = 0., d, INCR 
DISTANCE = X/d 
IF (DISTANCE.GE.0.0.AND.DISTANCE.LE.0.5) THEN 
FREAL = FRSTREAL(X) 
FIMAG = FRSTIMAG(X) 
ELSE 
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FREAL = SCNDREAL(X) 
FIMAG = SCNDIMAG(X) 
ENDIF 


ZET = CMPLX(FREAL, FIMAG) 
AMAGNIT = CABS(ZET) 
PHASE = PHASER(FREAL, FIMAG) 
WRITE(14,20) DISTANCE, FREAL, FIMAG, AMAGNIT, PHASE 
FORMAT(1X, F12.6, 4X, F12.6, 4X, F12.6, 4X, F12.6, 
4X, F12.6) 

CONTINUE 

STOP 

END 


RaKkKKKKKKKKKKKKKKKKKKKKKE 


REAL FUNCTION FRSTREAL(X) 
REKRAEKEERERKRKERKKKRKRKKKKRK 


REAL K, KC3 
COMPLEX J, POWER, FUSCAT, FUINC, UI, FUNC 
COMMON K, H, a, THETAO 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
A = - 4.%H/d 

C = SQRT(1.+ A*A) 
Z = H + AtX 

J = CMPLX(0.,1.) 


DI = (STHO)**2.+ A*CTHO*(2.*STHO - A*CTHO) 
FUSCAT = J#(1. - DI/(2.*C#C)) 
FUINC = J*((A*STHO + CTHO)/C) 


POWER = - J*K*(X*STHO - Z*CTHO) 
UI = CEXP(FOWER) 

FUNC = - J*UI*(FUSCAT + FUINC) 
FRSTREAL = REAL(FUNC) 

RETURN 

END 


RRREKKKKKKKKKKKRKKKKKKEKKRKK 


REAL FUNCTION FRSTIMAG (X) 
RRAREKRKEKEKRKKEEKEKREREREREREE 


REAL K, KC3 
COMPLEX J, POWER, FUSCAT, FUINC, UI, FUNC 
COMMON K, H, a, THETAO 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 

A = - 4.*H/d 

C = SQRT(1.+ A*A) 

Z =H + A*X 

J = CMPLX(0.,1.) 

DI = (STHO)**2.+ A*CTHO*(2.*STHO - A*CTHO) 
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FUSCAT = J*(1. - DI/(2.*C#C)) 
FUINC = J*((A*STHO + CTHO)/C) 
POWER = - J*#K*(X*STHO - Z*CTHO) 
UI = CEXP(POWER) 

FUNC = - J*#UI*(FUSCAT + FUINC) 
FRSTIMAG = AIMAG(FUNC) 

RETURN 

END 


REKRKKKKKRKKKKKKKKKEKKKKKEK 


REAL FUNCTION SCNDREAL(X) 
REAKKERKREKKEEEKRERKKKRKKKKE 


REAL K, KC3 
COMPLEX J, POWER, FUSCAT, FUINC, UI, FUNC 
COMMON K, H, d, THETAO 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
A = + 4.*H/d 

C = SQRT(1.+ A*A) 
Z = -3*H + A*X 

J = CMPLX(0.,1.) 


DI = (STHO) **2.+ A*CTHO*(2.*STHO - A*CTHO) 
FUSCAT = J*(1. - DI/(2.*C#C) ) 
FUINC = J*((A*STHO + CTHO)/C) 


POWER = - J*K*(X*STHO - Z*CTHO) 
UI = CEXP(POWER) 

FUNC = - J*UI*(FUSCAT + FUINC) 
SCNDREAL = REAL (FUNC) 

RETURN 

END 


RkRKKKKKKRKKKKKKKKKKKKKKKEK 


REAL FUNCTION SCNDIMAG (X) 
RREREKREREREREEEKEEEEREREEK 


REAL K, KC3 
COMPLEX J, POWER, FUSCAT, FUINC, UI, FUNC 
COMMON K, H, d, THETAO 


STHO 
CTHO 
A 
Cc 


SIN (THETAO) 

COS (THETAO) 

+ 4.*H/d 

SQRT(1.+ A*A) 

Z -3*H + A*X 

J CMPLX(0.,1.) 

DI = (STHO)**2.+ A*CTHO*(2.*STHO - A*CTHO) 
FUSCAT = J*(1. - DI/(2.*C#C) ) 
FUINC = J*((A*STHO + CTHO)/C) 
POWER = - J*K*(X*STHO - Z*CTHO) 
UI = CEXP( POWER) 
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FUNC = - J*UI*(FUSCAT + FUINC) 
SCNDIMAG = AIMAG(FUNC) 

RETURN 

END 
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PROGRAM STRIPMON 


THIS PROGRAM CALCULATES THE MONOSTATIC ECHO WIDTH 
OF A THIN, PLANAL, CONDUCTING STRIP OF LENGTH 2*L. 


NPUT DATA: 

F = FREQUENCY IN GHz 

KOL = KO*L IS THE PRODUCT OF THE STRIP HALF 
LENGTH (L) TIMES KO (KO=2*PI/WAVELENGTH) 

INCLIN = STRIP INCLINATION (IN DEGREES) WITH RESPECT 
TO THE HORIZONTAL (+X) AXIS 

STEP = ANGLE OF INCIDENSE INCREMENT 

OUTPUT DATA: 


ECHOMONO = NORMALIZED - WITH RESPECT TO THE WAVELENGTH - 
MONOSTATIC ECHO WIDTH AS A FUNCTION OF THE 

INCIDENT ANGLE THETAO, WHICH HERE (MONOSTATIC 

CASE) IS EQUAL TO: 2*PI - ANGLE OF INCIDENCE. 


REAL KOL, INCLIN 

OPEN (UNIT=1, FILE='STRIPMON.DAT', FORM='FORMATTED', 
STATUS=' UNKNOWN ' ) 

READ (5, *) F,KOL, INCLIN, STEP 


PI = 4.*ATAN(1.) 
RAD = PI/180. 

A = TAN(INCLIN) 
C2 = 1. + A*A 


DO 100 QI = 0., 90., STEP 
THETAO = QI*RAD 
THETAR = 2.*PI - THETAO 


STHO = SIN(THETAO) 
CTHO = COS (THETAO) 
STH = SIN(THETAR) 
CTH = COS(THETAR) 


D = (STHO)**2. + A*CTHO*(2.*STHO - A*CTHO) 
= STH - STHO + ALPHA*(CTH + CTHO) 


ECHOMO=C2/(2.*PI) *(2.*(1.-D/(2.*C2) ) *SIN(KOL*G) /G) **2 
WRITE(1,*) QI,ECHOMO 

CONTINUE 

STCP 

END 


66 





qagagaagqaanqanaanaaagqagagaaaaaana 


100 





PROGRAM SsTRIPBI 


THIS PROGRAM CALCULATES THE BISTATIC ECHO WIDTH 
OF A THIN, PLANAL, CONDUCTING STRIP OF LENGTH 2*L. 


INPUT. DATA: 


F FREQUENCY IN GHz 


KOL KO*L IS THE PRODUCT OF THE STRIP HALF 
LENGTH (L) TIMES KO (K0=2*PI/WAVELENGTH) 
INCLIN = STRIP INCLINATION (IN DEGREES) WITH RESPECT 


TO THE HORIZONTAL (+X) AXIS 

INCID = ELECTRIC FIELD ANGLE OF INCIDENCE (IN DEGREES) 
WITH RESPECT TO THE VERTICAL (+2) AXIS 

STEP = OBSERVATION ANGLE INCREMENT 


OUTPUT DATA: 


ECHOBI = NORMALIZED BISTATIC ECHO WIDTH AS A 
FUNCTION OF THE OBSERVATION ANGLE, THETA, 
FOR A GIVEN INCIDENT ANGLE 


REAL KOL, INCLIN, INCID 

OPEN (UNIT=2, FILE='STRIPBI.DAT', FORM='FORMATTED', 
STATUS=' UNKNOWN ' ) 

READ(5,*) F,KOL,INCLIN, INCID,STEP 


PI = 4.*ATAN(1.) 
RAD = PI/180. 
THETAO = INCID*RAD 
TAN (INCLIN) 
1. + A*A 


DO 100 QI = -180, 180, STEP 
THETAD = QI 
THETAR = QI*RAD 

SIN (THETAO) 

COS (THETAO) 

SIN ( THETAR) 

COS (THETAR) 


Q 
4 
ro) 
(o) 
iow it ol 


D = (STHO)**2. + A*CTHO*(2.*STHO - A*CTHO) 
= STH - STHO + ALPHA*(CTH + CTHO) 


ECHOBI=C2/ (2. *PI) *(2.*(1.-D/(2.*C2) ) *SIN(KOL*G) /G) **2 
WRITE(2,*) THETAD, ECHOBI 
CONTINUE 
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PROGRAM SINUSMON 


THIS PROGRAM CALCULATES THE MONOSTATIC ECHO WIDTH OF A 
SINUSOIDAL SURFACE z(x) = H*cos([2*pi/d]*x) 

INPUT. DATA: 

FREQ = FREQUENCY IN GHz 

H, @ = PARAMETERS OF THE SINUSOIDAL SURFACE 

KOL = LIMIT OF INTEGRATION (KOL = K*L), WHERE : 


"K" IS THE PROPAGATION CONSTANT AND 
"EL" IS THE ON THE X-AXIS SURFACE PROJECTION 


LENGTH 
STEP = ANGLE OF INCIDENSE INCREMENT 
M = NUMBER OF EQUIDISTANT POINTS BETWEEN THE TWO 


NORMALIZED LIMITS (0 AND 1) OF INTEGRATION 
(M ODD, M < 999). 


OUTPUT DATA: 


ECHOMO = NORMALIZED (WITH RESPECT TO THE WAVELENGTH) 
MONOSTATIC ECHO WIDTH, AS A FUNCTION OF THE 
INCIDENT ANGLE (THETAO) WHICH HERE (MONOSTATIC 
CASE) IS EQUAL TO: 2*PI - OBSERVATION ANGLE. 


ECHO = NORMALIZED (WITH RESPECT TO THE SURFACE LENGTH 
2*L) MONOSTATIC ECHO WIDTH AS A FUNCTION OF THE 
ANGLE OF INCIDENCE, THETAO. 


REAL F(999) 

REAL K,KOL, KOLU, INCR, INCID,MAGN2_R 

COMPLEX J, R 

EXTERNAL FILON, FREAL, FIMAG 

COMMON K, H, W, THETAO, THETAR 

OPEN (UNIT=7, FILE='SINUSMON.DAT', FORM='"FORMATTED', 
STATUS=' UNKNOWN ' ) 

READ(5,*) FREQ, H, da, KOL, STEP, M 


PI = 4.*ATAN(1.) 
W = 2.*PI/d 

RAD = PI/180. 

K = 20.*PI*FREQ/3. 
INCR = 1./M 
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DO 100 QI = 0., 90., STEP 
THETAO = QI*RAD 
THETAR = 2.*PI - THETAO 
TSO = KOL*(SIN(THETAR) - SIN(THETAO) ) 
SGN = SIGN(1,TSO) 
ABSTSO = ABS(TSO) 


Cc 
I=il1 
DO 40 U = 0., 1., INCR 
KOLU = KOL*U 
F(I) = MOEVREAL(KOLU) 
I=Ie+i 
40 CONTINUE 
Rl = FILON(F,ABSTSO,0.,1.,M,999,1) 
Cc 
I=1l 
Do 45 U = 0., 1., INCR 
KOLU = KOL*U 
F(I) = MOEVIMAG(KOLU) 
I=I+i1 
45 CONTINUE 
R2 = FILON(F,ABSTSO,0.,1.,M,999,1) 
Cc 
I=il 
DO 50 U = 0O., 1., INCR 
KOLU = KOL*U 
F(I) = MOODREAL(KOLU) 
I=I+i1 
50 CONTINUE 
R3 = FILON(F,ABSTSO,0.,1.,M,999,2) 
Cc 
I=il 
DO 55 U = 0., 1., INCR 
KOLU = KOL*U 
F(I) = MOODIMAG(KOLU) 
Ir=Ir+e+i1 
55 CONTINUE 
R4 = FILON(F,ABSTSO,0.,1.,M,999,2) 
Cc R : calculated value of the integral 
Cc 
R = (R1 - SGN*R4) + J*(R2 + SGN*R3) 
RR = Ri - SGN*R4 
RI = R2 + SGN*R3 
MAGN2_R = RR*RR + RI*RI 
= ECHO = sigma/(2*L) 
ECHO = MAGN2_R * KOL/8. 
Cc ECHOMO = sigma/wavelength 
ECHOMO = ECHO*KOL/PI 
Cc 


WRITE(7,*) QI, ECHOMO, ECHO 
100 CONTINUE 
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STOP 
END 


REKKEKKREKEEKREKKEREEKERKKKKKKKKKKRKKKKKKKKKKKKKS 


REAL FUNCTION FILON (F, T, A, B, M, MAX, KEY) 
HRERERKERRREKEEEERERKEREEREEEEREREEEEREEREREREREKEE 


THIS PROGRAM COMPUTES SP LF (x) *cos(Tx) Jax IF "KEY" = 1, 
AND f, (F(x) *sin(Tx)]dx WHEN "KEY" = 2, USING FILON'S 
METHOD. FUNCTION "F(x)" IS ASSUMED TO BE GIVEN AS A 
TABLE OF FUNCTIONAL VALUES IN THE ARRAY AT "M" 
EQUIDISTANT POINTS FROM "A" To “Bt, "M" ODD. 

"MAX" IS THE DIMENSION OF THE ARRAY DECLARED IN THE 
CALLING PROGRAM (FROM [REF. 9]). 


INTEGER KEY, M, N, I, MAX 

REAL F(MAX), AL, BE, GA, SUM, TH, H, S, C, S1, Cl, S2, C2 
REAL Fl, F2, SU, SU1, Al, A, B, T, TH2, TH4 

N=M-1 

H = (B - A) / FLOAT (N) 

TH = T * H 

S = SIN (TH) 


IF (ABS (TH) .LE. 0.1) THEN 

TH4 = TH2 * TH2 

= 2. / 45. * TH * TH2 * (1.- TH2 / 7. + TH4 / 105.) 
BE = 2. / 3. + 2. / 15. * TH2 - 4. / 105. * TH4 


GA = 4. / 3. - 2. / 15. * TH2 + 1. / 210. * TH4 
ELSE 
AL=1. / TH +S * C / TH2 - 2. * S * S / (TH2 * TH) 


= Ss 
BE = (2. +2. * C* C - 4. %* 5S %*C/ TH) / TH2 


GA = 4. * (S / TH - C) / TH2 
END IF 
SUM = 0.0 
Fl = F (1) 
F2 = F (M) 
S1 = SIN (A * T) 
S2 = SIN (B * T) 
Cl = COS (A * T) 
C2 = COS (B * T) 
Al=A+H 


SU1 = -0.5 * (F2 * C2 - Fl * Cl) 
2 2 


= td N, 
SUM = SUM + F (I) * COS (Al * T) 
Al = Al + H 
SU1 = SU1 + F (I+1) * COS (Al * T) 
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GO TO 4 

IF (T .EQ. 0.) THEN 
FILON = 0. 
RETURN 

END IF 


SU = Fl * Cl - F2 * C2 
SU1 = -0.5 * (F2 * S2 - Fl * S1) 
DO 5 I = 2, N, 2 
SUM = SUM + F (I) * SIN (Al * T) 
Al = Al +H 
SU1 = SU1 + F (I+1) * SIN (Al * T) 
Al = Al +H 
FILON = H * (AL * SU + BE * SU1 + GA * SUM) 


RAEKKKEKEREKKEEKEKEKKKKKKKKEE 


REAL FUNCTION MOEVREAL(RO) 
RRERKKREKKREREEKERERKKKKEKKEEE 


REAL K 
COMPLEX J, BJC3, DJE, FE, EV, FACTOR 
COMMON K, H, W, THETAO, THETA 


SIN (THETAO) 
COS (THETAO) 
SIN (THETA) 
COS (THETA) 
W*RO/K 

SIN (WROK) 

= COS (WROK) 

J = CMPLX(0.,1.) 
A = -H*W*Sw 
B= 
B 


-H*W*W/K*Cw 
A= ABS (B) 
C = SQRT(1.+ A*A) 
D = (STHO)*(STHO) - (A*CTHO) * (A*CTHO) 
E = 


B*CTHO/ (C#C) 
C3 = C*CKC 
BJC3 = BA + J*C3 
DIE = D + J#E 


TSE = K*H*CwW* (CTH - CTHO) 
FACTOR = 2.*C®CEXP(J*TSE) 


FE 


BA/ (2.*C3) *(1-BA/(4.*BJC3)) + C/(2.*BJC3) *DJE 
+J* (1+CTHO/C) 

EV = FACTOR*FE 

FEVREAL = REAL(EV) 
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REAL FUNCTION MOEVIMAG (RO) 
akhkkkkhkkeeereeeeeeerkkhret 


REAL K 
COMPLEX J, BJC3, DJE, FE, EV, FACTOR 
COMMON K, H, W, THETAO, THETA 


SIN (THETAO) 
COS (THETAO) 
SIN (THETA) 
COS (THETA) 
W*RO/K 
SIN (WROK) 
= COS (WROK) 
CMPLX(0.,1.) 
-H*W*SW 
~H*W*W/K*Cw 
ABS (B) 
SQRT(1.+ A*A) 
(STHO)*(STHO) - (A*CTHO) *(A*CTHO) 
B*CTHO/ (C*C) 
= CRCKC 
= BA + J*C3 
DIE = D+ J*E 


Cae g ame 
WwW 
Huuiuda 


TSE = K*H*CW*(CTH - CTHO) 
FACTOR = 2.*C*CEXP(J*TSE) 


FE = BA/(2.*C3)*(1-BA/(4.*BJC3)) + C/(2.*BJC3) *DJE 
+J* (1+CTHO/C) 

EV = FACTOR*FE 

FEVIMAG = AIMAG(EV) 

RETURN 

END 


RERKRKKKKKKKKKRKKKKKKKKKKKKK 


REAL FUNCTION MOODREAL(RO) 
HRRRRAKEREREEEREKREEKEKKKEKEE 


REAL K 
COMPLEX J, BJC3, DJE, FO, OD, FACTOR 
COMMON K, H, W, THETAO, THETA 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
STH = SIN(THETA) 

CTH = COS(THETA) 
WROK = W*RO/K 

SW = _SIN(WROK) 
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= COS (WROK) 
CMPLX(0.,1.) 
-H*W*SW 
-H*WAW/K*CwW 
ABS (B) 
SQRT(1.+ A*A) 
2. *A*STHO*CTHO 
A*B*STHO/ (C*C) 
C3 = CeCEC 

BJC3 = BA + J*C3 
DIE =D + J*E 


HOO PAG 


TSE = K*H*CW* (CTH - CTHO) 
FACTOR = 2.*C*CEXP(J*TSE) 


FO = C/(2.*BJC3)*DJE + J*A*STHO/C 
OD = FACTOR*FO 

FODREAL = REAL(OD) 

RETURN 

END 


RREKKEREKKREKKKKKKKKKKKKEEK 


REAL FUNCTION MOODIMAG (RO) 
KARKKRARRKREKEKEKRERKEKEEKK 


REAL K 
COMPLEX J, BJC3, DJE, FO, OD, FACTOR 
COMMON K, H, W, THETAO, THETA 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
STH = SIN(THETA) 
CTH = COS(THETA) 
WROK = W*RO/K 

SW = SIN(WROK) 
CW = COS(WROK) 

J = CMPLX(0.,1.) 

A = -H*W*SW 

B = -H*W*W/K*CW 
BA= ABS(B) 

C = SQRT(1.+ A*A) 
D = 2.*A*STHO*CTHO 
E = A*B*CTHO/(C*C) 


C3) = C®CRC 


TSE = K*H*CW*(CTH - CTHO) 
FACTOR = 2.*C*CEXP(J*TSE) 


FO = C/(2.*BJC3)*DJE + J*A*STHO/C 
OD = FACTOR*FO 
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FODIMAG = AIMAG(OD) 
RETURN 
END 


qgqgaggqNgNQqNagNNANNNNANANNNNNNNANANNANNNNNNNNNANAN 


PROGRAM SINUSBI 


THIS PROGRAM CALCULATES THE BISTATIC ECHO WIDTH OF A 
SINUSOIDAL SURFACE z(x) = H*cos([{2*pi/d]*x) 


NPU A: 

FREQ = FREQUENCY IN GHz 

INCID = ELECTRIC FIELD ANGLE OF INCIDENSE (IN DEGREES) 
WITH RESPECT TO THE VERTICAL (+Z) AXIS 

H, a@ = PARAMETERS OF THE SINUSOIDAL SURFACE 

KOL = LIMIT OF INTEGRATION (KOL = K*L), WHERE : 


"K" IS THE PROPAGATION CONSTANT AND 
"L" IS THE ON THE X-AXIS SURFACE PROJECTION 
LENGTH 


STEP = OBSERVATION ANGLE INCREMENT 


M = NUMBER OF EQUIDISTANT POINTS BETWEEN THE TWO 
NORMALIZED LIMITS (0 AND 1) OF INTEGRATION 
(M ODD, M < 999). 


OUTPUT DATA: 


ECHOBI = NORMALIZED (WITH RESPECT TO THE WAVELENGTH) 
BISTATIC ECHO WIDTH AS A FUNCTION OF THE 
OBSERVATION ANGLE THETA. 


ECHO = NORMALIZED (WITH RESPECT TO THE SURFACE LENGTH 
2*L) BISTATIC ECHO WIDTH AS A FUNCTION OF THE 
OBSERVATION ANGLE THETA. 


REAL F(999) 

REAL K,KOL, KOLU, INCR, INCID,MAGN2_R 

COMPLEX J, R 

EXTERNAL FILON, FREAL, FIMAG 

COMMON x, H, W, THETAO, THETAR 

OPEN (UNIT=8, FILE='SINUSBI.DAT', FORM='FORMATTED', 
STATUS=' UNKNOWN ' ) 


READ(5,*) FREQ, INCID, H, d, KOL, STEP, M 
PI = 4.*ATAN(1.) 
W = 2.*PI/d 


RAD = PI/180. 
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40 


45 


50 


55 


THETAO = INCID*RAD 
K = 20.*PI*FREQ/3. 
INCR = 1./M 


= -90., 90., STEP 
THETAD = QI 

THETAR = THETAD*RAD 

TSO = KOL*(SIN(THETAR) - SIN(THETAO) ) 
SGN = SIGN(1,TSO) 

ABSTSO = ABS(TSO) 


I=1 
DO 40 U = 0., 1., INCR 
KOLU = KOL*U 
F(I) = BIEVREAL(KOLU) 
I=I+1 
CONTINUE 
Rl = FILON(F,ABSTSO,0.,1.,M,999,1) 
I=1 
DO 45 U = 0., 1., INCR 
KOLU = KOL*U 
F(I) = BIEVIMAG(KOLU) 
I=I+i 
CONTINUE 
R2 = FILON(F,ABSTSO,0.,1.,M,999,1) 
I=1 
DO 50 U = 0., 1., INCR 
KOLU = KOL*U 
F(I) = BIODREAL(KOLU) 
I=nt+i 
CONTINUE 
R3 = FILON(F,ABSTSO,0.,1.,M,999,2) 
I=1 
DO 55 U = 0., 1., INCR 
KOLU = KOL*U 
F(I) = BIODIMAG(KOLU) 
I=I+1 
CONTINUE 


R4 = FILON(F,ABSTSO,0.,1.,M,999,2) 
R : calculated value of the integral 


R = (Rl - SGN*R4) + J*(R2 + SGN*R3) 
RR = R1 - SGN*R4 
RI = R2 + SGN*R3 


MAGN2_R = RR*RR + RI*RI 

ECHO = sigma/(2*L) 
ECHO = MAGN2_R * KOL/8. 

ECHOBI = sigma/wavelength 
ECHOBI = ECHO*KOL/PI 


77 


WRITE(8,*) QI, ECHOBI, ECHO 
100 CONTINUE 


STOP 
END 
Cc REKRKRKRKRKKKKKKKKKKKKKKRKKRKKK 
REAL FUNCTION BIEVREAL(RO) 
Cc REKRKKEKKEKKRKKKKKRKKRKKKRKRKKKKK 
Cc 
REAL K 
COMPLEX J, BJC3, DJE, FE, EV, FACTOR 
COMMON K, H, W, THETAO, THETA 
Cc 
STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
STH = SIN(THETA) 
CTH = COS(THETA) 
WROK = W*RO/K 
SW = SIN(WROK) 
CW = COS(WROK) 
J = CMPLX(0.,1.) 
A = -H*W*SW 
B = -H*W*W/K*CW 
BA= ABS (B) 
C = SQRT(1.+ A*A) 
D = (STHO)*(STHO) - (A*CTHO) *(A*CTHO) 
E = B*CTHO/ (CC) 
C3 = C*CKC 
BJC3 = BA + J*C3 
DIE = D+ J*E 
Cc 
TSE = K*H*CW*(CTH - CTHO) 
Cc 
FACTOR = 2.*C*CEXP(J*TSE) 
Cc 
FE = BA/(2.*C3) *(1-BA/(4.*BJC3)) + C/(2.*BJC3) *DJE 
& +J* (14+CTHO/C) 
EV = FACTOR*FE 
Cc 
FEVREAL = REAL(EV) 
RETURN 
END 
Cc 
Cc KRKRAKRKKKRKKKKRKKARKRKRKRKKRREKRE 
REAL FUNCTION BIEVIMAG (RO) 
Cc RAREKRKRKKRKKKKKRKKRKRKRKKKRKRKRKRREE 
Cc 
REAL K 
COMPLEX J, BJC3, DJE, FE, EV, FACTOR 
COMMON K, H, W, THETAO, THETA 
Cc 


STHO = SIN(THETAO) 
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CTHO = COS(THETAO) 
STH = SIN(THETA) 
CTH = COS(THETA) 
WROK = W*RO/K 

SW = SIN(WROK) 
CW = COS(WROK) 

J = CMPLX(0.,1.) 
A = -H*W*SW 

B = -H*W*W/K*CW 
BA= ABS (B) 

C = SQRT(1.+ A*A) 
D = (STHO)*(STHO) - (A*CTHO) *(A*CTHO) 
E = B*CTHO/(C#C) 
C3 = CeCKC 

BJC3 = BA + J*C3 
DIE = D + J*E 


TSE = K*H*CW*(CTH - CTHO) 
FACTOR = 2.*C*CEXP(J*TSE) 


FE = BA/(2.*C3) *(1-BA/(4.*BJC3)) + C/(2.*BJC3) *DJE 
+J* (1+CTHO/C) 

EV = FACTOR*FE 

FEVIMAG = AIMAG(EV) 

RETURN 

END 


REKRRAKRKERERKERKEKKKKERKEKKKKE 


REAL FUNCTION BIODREAL(RO) 
HRRKKKEEKKEEKREKKEKKKKKKKKKKKEE 


REAL K 
COMPLEX J, BJC3, DJE, FO, OD, FACTOR 
COMMON K, H, W, THETAO, THETA 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
STH = SIN(THETA) 
CTH = COS(THETA) 
WROK = W*RO/K 

SW = SIN(WROK) 
CW = COS(WROK) 

J = CMPLX(0.,1.) 

A = -H*W*SW 

B = -h*W*W/K*Cw 
BA= ABS (B) 

C = SQRT(1.+ A*A) 
D = 2. *A*STHO*CTHO 
E = -A*B*STHO/(C*C) 
C3 = C*CKC 

BJC3 = BA + J*C3 
DIE =D + JkE 
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TSE = K*H*CW*(CTH - CTHO) 
FACTOR = 2.*C*CEXP(J*TSE) 


FO = C/(2.*BJC3)*DJE + J*A*STHO/C 
OD = FACTOR*FO 

FODREAL = REAL(OD) 

RETURN 

END 

HEKKKEEEREEKEREREKKRKKKKKKEKE 


REAL FUNCTION BIODIMAG (RO) 
HRKKKERAREKREKREKEREKKEERKEKE 


REAL K 
COMPLEX J, BJC3, DJE, FO, OD, FACTOR 
COMMON K, H, W, THETAO, THETA 


STHO = SIN(THETAO) 
CTHO = COS(THETAO) 
STH = SIN(THETA) 
CTH = COS(THETA) 
WROK = W*RO/K 

SW = SIN(WROK) 
CW = COS(WROK) 

J = CMPLX(0.,1.) 

A = -H*W*SW 

B = -H*W*W/K*CW 
BA= ABS (B) 

C = SQRT(1.+ A*A) 
D = 2.*A*STHO*CTHO 
E = A*B*STHO/ (C*C) 
C3 = CCRC 

BJC3 = BA + J#C3 
DIE = D + J*E 


TSE = K*H*CwW*(CTH - CTHO) 
FACTOR = 2.*C*CEXP(J*TSE) 


FO C/(2.*BJC3) *DJE + J*A*STHO/C 
oD FACTOR* FO 

FODIMAG = AIMAG(OD) 

RETURN 
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